% =========================================================================
% run_simulation.m - (高效验证版) NMPC能量管理仿真主程序
%
% 核心修改:
% 1. 移除了waitbar进度条，让仿真以MATLAB的最快速度运行。
% 2. 增加了tic/toc计时功能，用于显示总仿真耗时。
% 3. 所有其他的仿真逻辑、参数、扰动生成方式均保持您原始版本不变。
% =========================================================================

%% 1. 初始化
InitParameters; % 运行初始化脚本，创建所有参数和nlobj控制器

%% 2. 设置初始条件和扰动序列 (保持您的原始设定)
x0 = [0.8; FCParams.P_fc_min; 0]; % 初始状态: [SOC=80%, P_fc=min, alpha=0]
u0 = 0.5;                         % 初始控制输入

% 创建一个模拟的负载功率序列 (P_load)
N_sim = floor(SimParams.T_sim / SimParams.Ts); % 总仿真步数
P_load_sequence = 150 * ones(N_sim, 1);
P_load_sequence(50:150) = 400; % 5-15s: 高功率需求
P_load_sequence(151:175) = linspace(400,200,25);
P_load_sequence(176:300) = 200;% 15-30s: 中等功率
P_load_sequence(301:310) = linspace(200,450,10);
P_load_sequence(310:400) = 450;% 31-40s: 再次高功率
P_load_sequence(401:410) = linspace(450,150,10);
P_load_sequence(411:470) = 150; % 40-50s: 低功率
P_load_sequence(471:500) = linspace(150,0,30);
% P_load需要加上噪声
noise_sequence = 1+0.2*randn(size(P_load_sequence));
P_load_sequence = P_load_sequence.*noise_sequence;

%% 3. 准备仿真历史记录 (保持您的原始设定)
x_history = zeros(N_sim+1, nx);
u_history = zeros(N_sim, nu);
x_history(1,:) = x0';
xk = x0; % 当前状态

%% 4. 仿真主循环 (修改运行框架)
disp('开始仿真... (程序将全速运行，请稍候)');

nloptions = nlmpcmoveopt; % 创建NMPC求解器选项
tic;
for k = 1:N_sim
    % --- a. 获取当前扰动 ---
    P_load_k = P_load_sequence(k);
    
    % --- b. 【核心】调用NMPC求解器计算最优控制量 ---
    nloptions.Parameters = {P_load_k};
    [uk, nloptions, info] = nlmpcmove(nlobj, xk, u_history(max(1,k-1),:), [], [], nloptions);
    
    % --- c. 将最优控制量应用到被控对象 ---
    xk_next = prediction_model(xk, uk, P_load_k);
    
    % --- d. 保存历史数据并更新状态 ---
    u_history(k,:) = uk';
    x_history(k+1,:) = xk_next';
    xk = xk_next;
    
    % 【修改】移除了此处的 waitbar(k/N_sim, hbar);
end

sim_duration = toc; % 【修改】结束计时

disp('仿真结束。');
disp(['总耗时: ', num2str(sim_duration), ' 秒。']); % 【修改】显示总耗时

%% 5. 结果可视化 (保持您的原始设定)
% 创建时间轴
t = 0:SimParams.Ts:SimParams.T_sim;

% 绘制状态量
figure('Name', '状态量变化曲线');
subplot(3,1,1);
plot(t, x_history(:,1), 'b', 'LineWidth', 2);
title('SOC 变化曲线'); grid on; ylabel('SOC');
ylim([0.15, 0.85]);

subplot(3,1,2);
plot(t, x_history(:,2), 'r', 'LineWidth', 2);
title('燃料电池输出功率'); grid on; ylabel('P_{fc} (W)');

subplot(3,1,3);
plot(t, x_history(:,3), 'g', 'LineWidth', 2);
title('燃料电池累积衰减'); grid on; ylabel('\alpha_{decay} (%)');
xlabel('时间 (s)');

% 绘制控制量和扰动
figure('Name', '控制与功率分配');
subplot(2,1,1);
plot(t(1:end-1), u_history, 'm', 'LineWidth', 2);
title('控制量 \beta 变化曲线'); grid on; ylabel('\beta');
ylim([-0.1, 1.1]);

subplot(2,1,2);
plot(t(1:end-1), P_load_sequence, 'k--', 'LineWidth', 2);
hold on;
P_fc_out = u_history .* P_load_sequence;
P_bat_out = P_load_sequence - P_fc_out;
area(t(1:end-1), [P_fc_out, P_bat_out]);
title('功率分配'); grid on; ylabel('功率 (W)');
legend('总负载 P_{load}', '燃料电池 P_{fc}', '锂电池 P_{bat}');
xlabel('时间 (s)');